library(foreign)

dorm <- read.dta("dorm results.dta")

dorm$outgroup.affect <- rowMeans(dorm[,c("mex1","afam1","muslim1")])
dorm$sophomore <- as.numeric(dorm$class == 2)
dorm$junior <- as.numeric(dorm$class == 3)
dorm$senior <- as.numeric(dorm$class == 4)

colMeans(subset(dorm[,c("age_01","female","knowroommate","freshman","sophomore","junior","senior","outgroup.affect", "lib","mod","cons","george1_01","auth01","hhi_rev_01")], dorm$hhi_rev_01 > median(dorm$hhi_rev_01)), na.rm = TRUE)

x <- subset(dorm[,c("age","female","knowroommate","freshman","sophomore","junior","senior","outgroup.affect", "lib","mod","cons","george1_01","auth01","hhi_rev_01")], dorm$hhi_rev_01 > median(dorm$hhi_rev_01))

Sx <- cov(x, use = "complete.obs")

mean(mahalanobis(x, colMeans(x, na.rm = TRUE), Sx), na.rm = TRUE)

x <- subset(dorm[,c("age","female","knowroommate","freshman","sophomore","junior","senior","outgroup.affect", "lib","mod","cons","george1_01","auth01","hhi_rev_01")], hhi_rev_01 < median(dorm$hhi_rev_01))

Sx <- cov(x, use = "complete.obs")

mean(mahalanobis(x, colMeans(x, na.rm = TRUE), Sx), na.rm = TRUE)


